%	lans_zmoment	- Compute Zernike moments of a 2-D image 
%
%	[A_nm,zmlist,cidx[,V_nm]]	= lans_zmoment(img[,n[,m]])
%
%	_____OUTPUTS____________________________________________________________
%	A_nm	Zernike moments					(row vector)
%	zmlist	keeps order and repetition information		(Mx2 matrix)
%		M = total # of moments
%	cidx	1-D indices of retained unit circle region	(vector)
%	V_nm	Zernike polynomials size of img			(col vectors)
%		to save memory, exclude this output term
%
%	_____INPUTS_____________________________________________________________
%	img	row by col image of reals			(matrix)
%	n	order						(vector)
%		starting from 0
%	m	repetition					(integer)
%
%	_____NOTES______________________________________________________________
%	- for demo, call function without parameters 	
%
%	lans_zmoment(img)
%		computes moment of order 0 repetition 0 
%	lans_zmoment(img,n)
%		computes all moments of order specified in VECTOR n
%	lans_zmoment(img,n,m)
%		computes moment of VECTOR order n repetition m
%		- |m|<=n and n-|m| even
%	- img MUST first be normalized w.r.t. scale and translation
%		i.e. centroid of image is at image center, and image
%		has geometrical moment m00 = fixed number/scale
%	- x,y follows raster scanning convention
%	- rotation invariant in the magnitude, i.e. |A_nm|
%	- to start from a particular order, specify a range in n, e.g. [4 5 6]
%	- n=0 is useless for centered data as it is simply the center of mass
%
%	_____SEE ALSO___________________________________________________________
%	lans_zmrecon	lans_zpoly	lans_fac	lans_centroid	lans_clipimg
%
%	(C) 2000.11.27 Kui-yu Chang
%	http://lans.ece.utexas.edu/~kuiyu

%	This program is free software; you can redistribute it and/or modify
%	it under the terms of the GNU General Public License as published by
%	the Free Software Foundation; either version 2 of the License, or
%	(at your option) any later version.
%
%	This program is distributed in the hope that it will be useful,
%	but WITHOUT ANY WARRANTY; without even the implied warranty of
%	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%	GNU General Public License for more details.
%
%	You should have received a copy of the GNU General Public License
%	along with this program; if not, write to the Free Software
%	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
%	or check
%			http://www.gnu.org/

function	[A_nm,zmlist,cidx,V_nm]	= lans_zmoment(img,n,m)

if nargin>0
%__________ REGULAR ____________________________________________________________
if nargin==1
	n	= 0;
end

d		= size(img);			% dimension
img		= double(img);

%----------	Normalize co-ordinates to [-sqrt2,sqrt2] (centroid assumed @ origin)
xstep		= 2/(d(1)-1);
ystep		= 2/(d(2)-1);
[x,y]		= meshgrid(-1:xstep:1,-1:ystep:1);


%!!!!uprava moje
[x_,y_,x0,y0] = hwozd_getMassCenter(hwozd_Image(img));
x = x + x0;
y = y + y0;

%----------	compute unit disc mask
circle1		= x.^2 + y.^2;
inside		= find(circle1<=1);
mask		= zeros(d);
mask(inside)	= ones(size(inside));

%----------	mask pixels lying inside/on unit circle
[cimg,cidx]	= lans_clipimg(img,mask);
z		= lans_clipimg(x+i*y,mask);
p		= abs(z);
theta		= angle(z);

%----------	Compute Zernike moments
c	= 1;
for order=1:length(n)
	n1	= n(order);
%	disp(sprintf('computing order %d',n1));
	if nargin<3
		m	= lans_zpossible(n1);		% repetition unspecified
	end
	
	for r=1:length(m)
		V_nmt		= lans_zpoly(n1,m(r),p,theta);
		zprod		= cimg.*conj(V_nmt);
		A_nm(c)		= (n1+1)*sum(sum(zprod))/pi;
		zmlist(c,1:2)	= [n1 m(r)];
		if nargout==4
			V_nm(:,c)	= V_nmt;
		end
		c		= c+1;
	end
end
%__________ REGULAR ends _______________________________________________________
else
%__________ DEMO _______________________________________________________________
clf;clc;
disp('running lans_zmoment.m in demo mode');

dim		= 16;
c		= round(dim/2);
l		= 4;
r		= 3;
gimg		= zeros(dim);
gimg(c-l:c+l,c-l:c+l)	= ones(2*l+1);
gimg		= gimg*255;

gimg(c-r:c+r,c-r:c+r)	= zeros(2*r+1);

%----------	plot scale/translation normalized image
clf;
iimg=lans_invariant(gimg,'-scale 0 -translation 0');
image(iimg);
colormap gray;
centroid= lans_centroid(double(iimg));
xcen	= centroid(1);
ycen	= centroid(2);
hold on;
plot(xcen,ycen,'r+');
axis square;
title('Original Image');
drawnow;

%----------	compute Zernike Moments up to an order
order				= 0:30;
disp('Computing Zernike moments of order 30');
[A_nm,zmlist,cidx,V_nm]		= lans_zmoment(iimg,order);


%----------	compute/plot reconstructed image
disp('Reconstructing image from Zernike moments (order 1, 2, ..., 30)');
[reimg]	= lans_zmrecon(A_nm,V_nm,zmlist,cidx,[dim dim],50,1);
%__________ DEMO ends __________________________________________________________
end
